Glass transition and random close packing above three dimensions 
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Motivated by a recently identified severe discrepancy between a static and a dynamic theory of 
glasses, we numerically investigate the behavior of dense hard spheres in spatial dimensions 3 to 
12. Our results are consistent with the static replica theory, but disagree with the dynamic mode- 
coupling theory, indicating that key ingredients of high-dimensional physics are missing from the 
latter. We also obtain numerical estimates of the random close packing density, which provides new 
insights into the mathematical problem of packing spheres in large dimension. 



Studying physical theories upon varying the dimen- 
sionahty of space d has illuminated a variety of prob- 
lems over the last century. Unifying fundamental forces 
through general relativity and string theories immedi- 
ately come to mind, but so does the treatment of critical 
phenomena through the renormalization group approach. 
Similar contemplations may now benefit another grand 
challenge of condensed matter, that of describing glass 
formation [l|. Including d within the set of control pa- 
rameters nicely complements two- and three-dimensional 
experiments and simulations by, first, surmounting some 
of the technical challenges encountered in low dimen- 
sions [31; and, second, by bringing the problem in contact 
with rich mathematical fields, such as coding theory [Sf. 
Obtaining a complete and consistent high-dimensional 
theory of amorphous states would also crucially pave the 
way for a more systematic understanding of finite dimen- 
sional effects, and ultimately provide a reliable theory of 
the glass transition, jamming, and related phenomena. 

The Random First Order Transition (RFOT) theory 
has emerged as one of the main contenders in the quest 
for a complete description of the glass transition. Its 
foundation was posed in the late eighties, when it was 
realized that a class of mean-field spin models share the 
phenomenology of glass formers [J]. These models are 
very abstract and do not obviously present the micro- 
scopic features of particle-based systems, but they helped 
turn RFOT theory into a quantitative method for study- 
ing glasses, thanks to the development of the dynamic 
mode-coupling theory (MCT) [5[ and the static replica 
theory (RT) [6|. The two quantitative theories of the 
glass transition are intimately related [7|, |8|, but make 
use of different approximations in order to obtain closed- 
form structures [Sj. Diverging physical descriptions rc- 
suh. The predictions of MCT 0, [13] and RT [ul in- 
deed strongly disagree in the limit of large d, with MCT 
suggesting that the glass transition happens at densities 
much larger than that of random close packing provided 
by RT. This situation is particularly worrying [12] . Be- 
cause RFOT theory is based on a mean-field description. 



it should be exact in large d. If MCT is really the dy- 
namic part of RFOT theory, then MCT and RT should 
lead to identical predictions, as proposed in [vf. Whether 
MCT is re ally a mean-field theory, however, has been 
questioned [l3| . In order to test the quality of the approx- 
imations and to clarify the connection with the original 
mean-field models, one ought to understand the source 
of the discrepancy, and to obtain a controlled limit of 
particle systems in which RFOT becomes exact. 

Independently from this issue, packing spheres in high 
dimensions is intimately related to several important 
mathematical problems, notably in the context of sig- 
nal digitalization and of error correcting codes |3| . It has 
been suggested that for large d disordered pac king s may 
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be more efficient than lattice-based versions 
though lattice geometries are strongly d dependent, the 
fluid structure and properties have been suggested to be 
much less sensitive to d, once d > 3 (see e.g. Ref. [15[). A 
general understanding of disordered packings in arbitrary 
d would thus also clarify the density scaling of amorphous 
packings and their potential mathematical role. 

In this letter, we compare simulations of hard spheres 
up to d = 12 with the predictions of RT and MCT. 
The results suggest that while RT offers a satisfying de- 
scription of dimensionality MCT fails at the task, which 
calls for a modified formulation of the dynamical theory 
of RFOT. The results also provide the "random close- 
packed" density in several d and give a hint of the scaling 
of this quantity for large d, which should allow compar- 
isons with the results of other theoretical treatments [l6| . 

Numerical simulations - We employ a modified 
Lubachcvsky-Stillingcr algorithm to density a low- 
density gas of N identical hard spheres of diameter a, en- 
closed in a periodic box of volume V , by growing the par- 
ticles at a constant rate 7 = fj, reported here in standard 



reduced units [17(, |18[ . Time evolution stops when the 
system reaches a high reduced pressure p = PP/ p = 10'^ 
measured by rescaling the mechanical pressure P by the 
number density p = N /V and the inverse temperature 
/3, which is thermostated to unity. The packing frac- 



tion is ip = pVd{<j/2), where Vd{R) is the volume of a 
d-dimensional baU of radius R. Our event-driven molec- 



ular dynamics scheme complements Ref. [17| 's earlier im- 
plementation of cubic blocking of space with spherical 
nearest-neighbor lists [19| . Because the volume of a ball 
inscribed in a cube tends to zero with growing d, re- 
markable efficiency gains are obtained from considering 
collisions with fewer neighbors. Up to d = 10, particles 
are grown at rates as low as 7 = 3 x 10^^, while rates 
of 10~^ and 10""^ are attained in d = 11 and 12, respec- 
tively. Systems with N = 8000 are simulated for d < 9 
and larger ones for d = 10-12 [2l|. These sizes ensure 
that even when the system is in its densest state the 
box edge remains larger than 2cr, which prevents a par- 
ticle from ever having two direct contacts with another 
one. There are strong reasons to believe that although 
relatively small these N nonetheless provide a reliable 
approximation of bulk behavior. First, with increasing 
d the box edge becomes less representative of the overall 
box size. The largest diagonals are yd larger and there 
are many more diagonals than edges. Second, by anal- 
ogy to spin systems, mean-field arguments indicate that 
for d ^ 00, a hypercube of side two is sufficient to cap- 
ture the full thermodynamic behavior. Even at the criti- 
cal point, finite-size corrections are proportional to 1/N^ , 
where the exponent 6 is model dependent (e.g. 1/2 at the 
ferromagnetic transition) , and do not directly involve the 
edge length L [22|. Similar results hold for dimensions 
greater than the upper critical dimension, where the ex- 
ponents coincide with the mean field ones. Third, the 
fluid structure is expected to become uniform at ever 
smaller distances with increasing d [17|, |20[. Nearest- 
neighbor ordering should thus mainly be influenced by 
particles in contact or nearly so, with the rest of the fluid 
acting as a continuum. Indeed, in the fluid phase, finite 
volume corrections are proportional to the pair correla- 
tion h(L), and at fixed L, h{L) goes to zero exponentially 
with d [20|. The validity of these rationalizations, which 
are consistent with the decorrelation property of high d 
sphere packings proposed in |14l | , are satisfactorily tested 
by simulations in d = 



211 



Numerical results - The compression results for d = 9 
shown in Fig. [T] are representative of the behavior ob- 
served for all d > 3, where crystallization is not observed 
on the time scale accessible to present computers |15l.ll7l|: 
in d = 3 crystallization is observed at small 7 leading to 
a drop in the pressure at intermediate density |17| . The 
system first follows the equilibrium fiuid equation of state 
(EOS) at low density and falls out of equilibrium at high 
density. Beyond this point, the pressure increases faster 
than in the equilibrium fluid and ultimately diverges at 
packing fraction (pj{'y). A Carnahan-Starling form 
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FIG. 1: Different compactions oi N = 8000 particles in d = 9. 
With growing ip, the pressure first evolves like the fluid EOS 
then like a free volume EOS. Extrapolated threshold glass and 
theoretical ideal glass lines illustrate the subsequent analysis 
(see text for details). 



captures well the pressure growth with (p in the fluid 
regime (Fig. [Ij , provided, for each d, that one fits Ad to 
the data from the slowest compression rate available [23J . 
Note that the coefficients A^^ are not identical to A^f 
adjusted to recover the correct third virial coefficient [2^ 
(see [2l|), but the values are quite close, and in any case 
this contribution vanishes with increasing d. 

In the high-density non-equilibrium regime, com- 
paction runs with different 7 follow separate branches 
along which the pressure evolution is dominated by the 
expulsion of free volume [25l [26| . Upon approaching jam- 
ming, the pressure is well approximated by 



Pfv(7,v) 



d^j(7)[l-/(7)] 
(^j (7) - (/? 



(2) 



(1) 



where both f{j) and 'p>j{"f) are extracted from fitting the 
simulation data for p > p^nin (see the table in [21|). Very 
close to jamming (p > 10^), 7(7) can be interpreted as 
the fraction of "rattlers" present, but in the regime where 
the fluid first becomes non-ergodic, caging heterogeneity 
results in a larger effective / [25[. We find that with 
decreasing 7, /{j) converges to values of order 10%, with 
only a weak d dependence (Fig.[T|). 

Data analysis - The numerical results qualitatively 
agree with the RFOT scenario (see Ref. [llj for details). 
According to the theory, the glassy states for moderately 
small 7 should converge as a power law 7" to a "thresh- 
old" glass that eventually jams at ipth- The dynamical 
transition density ip^ [27| separates the equilibrium fluid 
from this glass [28|, but at much slower compaction rates, 
7 < exp(— d), RFOT theory also predicts that activated 
events allow the system to remain in equilibrium up to 
higher densities. In this regime we expect a crossover to 
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FIG. 2: RFOT critical densities obtained from simulations, 
RT, and MCT. In the inset, the larger d interval shows the 
approach to the asymptotic d — >■ cxj limit, which happens for 
d f» 60 in both RT and MCT. 



a logarithmic dependence of the glass EOS on 7, e.g., of 
the form 1 / 1 log(7) | . As long as crystallization remains 



suppressed [l5l [17[, an extrapolation to 7 — >■ in this 



regime allows to identify a second order phase transition 
to an "ideal" glass, at a density (/Sk that corresponds to 
the Kauzmann point. This ideal glass then jams at the 
glass close-packed (GCP) density vgcp- Hence, from this 
mean-field treatment, stable amorphous packings should 
exist in the interval tp € [(/3th, y'GCp] jilt . 

Based on this description, it is clear that the acti- 
vated regime of RFOT is inaccessible in high dimensions. 
We therefore focus on the power-law regime and on the 
threshold glass. In the inset of Fig.[Tl the reported /(7) 
and (/5j(7) are both linear as functions of 7" with a ~ 0.5. 
The exponent is expected to weakly depend on dimen- 
sionality, especially at low d, but this value is within the 
numerically reasonable range for all systems studied. Be- 
cause the subsequent analysis is rather insensitive to the 
precise value of a, for simplicity it is kept constant. Ex- 
trapolating the results to 7 = gives the parameters Lpth 
and /th for the threshold glass free volume EOS reported 
in Fig. [T] Its intersection with the fluid EOS then pro- 
vides a numerical estimate for (/jj- These values are very 
close to previous numerical estimates in d = 3 [29| and 
d = 4 y. and the results for (/Jth agree with the results 
of Refs. [iTl |30|, which demonstrates the coherence of 
our analysis in low diixiensions. The numerical results 
reported in Table |T] and plotted in Fig. [2] further show 
that the dimensional evolution of ipth and ipd is smooth. 
Interestingly, dimensions where the crystal structure is 
singularly dense, such as d = 8 and d ~ 12 [3j, do not 
present any echo of that singularity, which illustrates the 
smooth d dependence of the fluid structure. 
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(RT-fCS) (RT+CS) (MCT+PY) 
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0.571 


0.651 


0.618 


0.684 


0.516 
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0.401 


0.467 


0.432 


0.487 


0.371 
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0.267 


0.319 


0.289 


0.331 


0.254 
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0.172 


0.209 


0.189 


0.219 


0.169 
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0.106 


0.133 


0.120 


0.140 


0.109 
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0.0648 


0.0821 


0.0738 


0.0875 


0.0699 


9 


0.0385 


0.0498 


0.0447 


0.0535 


0.0434 


10 


0.0226 


0.0297 


0.0266 


0.0319 


0.0264 


11 


0.0131 


0.0174 


0.0155 


0.0187 


0.0158 


12 


0.0075 


0.0101 


0.00891 


0.0108 


0.00964 
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TABLE I: Numerical values of the critical densities obtained 
from the simulations extrapolated at 7 = and the RFOT 
theories. The uncertainty on the numerical results is on the 
last reported digit [2ll |. 



Theory - The numerical results can also be quanti- 
tatively compared with the analytical estimates for the 
RFOT critical densities. Although a reliable calculation 
of (pd and y'th has not yet been obtained within RT, the 
latter gives results for (^k and (/Jgcp 



ll|, taking the 



Carnahan-Starling EOS with parameter A^'^ as input. 
As we mentioned above, the fitted values A^ are slightly 
different, but that difference only weakly perturbs the 
results. Because A^ is only available for d < 12, we 
use instead A^^ from Ref . [ll| . As expected from RFOT 
theory, ipK > ^d and (fiocp > Vth, and the values follow 
a similar trend with d (Table |l|. 

The MCT analysis follows the approach of Refs. [i[l^ 
using the Pcrcus-Yevick (PY) structure factor calculated 
iteratively with a numerical Hankel transformation of or- 
der d/2 — 1. Using instead the hypernetted chain (HNC) 
input for the structure factor does not strongly affect the 
results. For the dynamical and static theories to be con- 
sistent, the transition density (pucT predicted by MCT, 
reported in Table U and plotted in Fig. [21 should coincide 
with ip^. It is not the case. The MCT transition at (^mct 
increases too fast with d and around d = 13 will likely 
cross the numerical estimate of tpth- This situation is 
paradoxical, because the liquid should fall out of equilib- 
rium well before jamming occurs. The usual suggestion 
that activated events may improve the consistency by in- 
creasing (fMCT would here make things worse. 

The inset in Fig. [2] presents (pucT and pK for larger di- 
mensions. Both curves approach the results obtained by 
neglecting all structure for the liquid [20[ , which amounts 
to using the van der Waals' expression p = 1 + 2'*^^(^ for 
the fluid EOS and the Mayer function/(r) = -e{a - r) 
for the direct correlation function [9|, llOf ■ This treatment 
is exact for d — > 00 , a fact that is also consistent with our 
simulation results for the fluid EOS. In this asyniptotic 
large d regime, MCT predicts pmct ^ 0.22d'^ 2''^ U\M^ 



while RT predicts ipK ~ dlog{d) 2 "^ [n|. The simulation 



results at intermediate d thus give supporting evidence 
that the scaling predicted by MCT is incorrect. 

Conclusions - Our numerical results exte nd p revious 
estimates of the glass transition density [2, l29j and of 
the amorphous packing density [17| up to d = 12. These 
results are obtained partly thanks to methodological im- 
provements jl9[ over the approach of Ref. 17 1. Reas- 



suringly, our procedure for the data analysis, which is 



grounded in the RFOT scenario [ll|, produces results 
that are consistent with estimates obtained through al- 
ternate routes. Two key features arise from the work. 
First, our numerical estimates of (pth are slightly smaller 
than the values of V'gcp predicted by RT and follow a 
similar trend, indicating that RT provides a reliable pre- 
diction for the jamming density of hard spheres in high 
dimensions. In particular, the random close-packed state 
of frictionless hard spheres, which is expected to be proto- 
col dependent [ll|,l3l|j should fall within the range of den- 
sities predicted by the RFOT scenario ip € [</5th, <y5GCp]- 
We thus expect our numerical results for iy9th to be good 
estimates of that density, irrespective of the protocol 
used. Second, our results for the dynamical (glass) tran- 
sition ipd are slightly smaller than the ideal glass (Kauz- 
mann) transition ipK predicted by RT, as it should be. 
But at the same time, the MCT prediction (pucT does 
not coincide with ip^, unlike what the RFOT scenario 
suggests. The two quantities additionally show very dif- 
ferent trends with dimension and the discrepancy in- 
creases at large d, which suggests that the standard MCT 
formulation [5|, |9|, |lO[ is not exact at large d. Echoing 
and amplifying the call of Ref. [10| , a new theory should 
be constructed, in order to interpolate between MCT in 
d = 3 and the large d mean-field regime. 

Whether this theory will be a minor modification of 
MCT or a completely new theory remains an open ques- 
tion. Indeed, mean-field models displaying a RFOT are 
described by MCT- like equations [i,!!], and MCT itself 
has been described as a typ e of Landau theory of the 
dynamical transition 



12 



. It is thus possible that 
a modified MCT with a similar structure might be ex- 
act in large d. Steps in this direction have been already 



taken (28|, |33|. Constructing a complete RFOT theory 
that includes both RT and a modified MCT would be an 
important intellectual advance, as it would be the first 
complete theory of the glass and jamming transitions for 
particle-based systems. This advance would open the 
way towards a more systematic construction of RFOT, 
through a large d expansion or a renormalization group 
approach. Some of the approximations underlying the 
construction of MCT might also be better justified, lead- 
ing to a fully predictive glass theory for d < 3. 

We thank J. Kurchan, R. Mari, K. Miyazaki and 
R. Schilling for stimulating discussions. PC acknowledges 
NSF support No. DMR-1055586. 



[4; 
[s; 

[6; 
ir 

[8 
[9 

[lo; 
[11 

[12 

[13 

[14 
[15 
[16 
[17 
[18 

[19 



[20 
[21 

[22' 
[23 



[24; 
[25; 
[26; 

[27; 



[28' 
[29; 



R W. Anderson, Science 267, 1615 (1995). 

P. Charbonneau, A. Ikeda, J. A. van Meel, and 

K. Miyazaki, Phys. Rev. E 81, 040501(R) (2010). 

J. H. Conway and N. J. A. Sloane, Sphere Packings, 

Lattices and Groups ( Springer- Verlag, New York, 1988); 

H. Cohn, Proceedings of the international congress of 

mathematicians (ICM 2010), 2416 (|arXiv: 1003.3053p . 

T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, 

Phys. Rev. A 40, 1045 (1989). 

W. Gotze, Complex Dynamics of Glass-Forming Liquids, 

International Series of Monographs on Physics, Vol. 143 

(Oxford University Press, Oxford, 2009). 

M. Mezard and G. Parisi, Phys. Rev. Lett. 82, 747 

(1999). 

T. R. Kirkpatrick and P. G. Wolynes, Phys. Rev. A 35, 

3072 (1987). 

G. Szamel, Europhys. Lett. 91, 56004 (2010). 

B. Schmid and R. Schilling, Phys. Rev. E 81, 041502 
(2010); R. Schilling and B. Schmid, Phys. Rev. Lett. 106, 
049601 (2011). 

A. Ikeda and K. Miyazaki, Phys. Rev. Lett. 104, 255704 
(2010); Phys. Rev. Lett. 106, 049602 (2011). 
G. Parisi and F. Zamponi, Rev. Mod. Phys. 82, 789 
(2010). 

J.-P. Bouchaud, "The mode-coupling theory of su- 
percooled liquids: Does it wear any clothes?" 
http://www.condmatjournalclub.org/?p=1022 (2010) . 
W. Kob and R. Schilling, J. Phys.: Condens. Matter 3, 
9195 (1991). 

S. Torquato and F. H. Stillinger, Exp. Math. 15, 307 
(2006). 

J. A. van Meel, B. Charbonneau, A. Fortini, and P. Char- 
bonneau, Phys. Rev. E 80, 061110 (2009). 

C. Song, P. Wang, and H. A. Makse, Nature 453, 629 
(2008); Y. Jin et al. Phys. Rev. E 82, 051126 (2010). 
M. Skoge, A. Donev, F. H. Stillinger, and S. Torquato, 
Phys. Rev. E 74, 041127 (2006). 

In contrast to the original implementation, a canonical 
definition of temperature is here used, which rescales the 
time unit ^J^rruj^ for particles of unit mass m. 
L. Verlet, Phys. Rev. 159, 98 (1967); A. Donev, 
S. Torquato, F. H. Stillinger, J. Comp. Phys. 202, 765 
(2005). 

G. Parisi and F. Slanina, Phys. Rev. E 62, 6554 (2000). 
See Appendix to this paper for numerical details, finite- 
size studies, and error analysis. 

G. Parisi, Statistical Field Theory (Perseus Group, 1998). 
The many virial coefficients known for hard spheres, e.g., 
Ref. [SJI , insufficiently capture the equation of state near 
the onset density of the non-ergodic regime. 
Y. Song, E. A. Mason, and R. M. Stratt, J. Phys. Chem. 
93, 6916 (1989). 

A. Donev, S. Torquato, and F. H. Stillinger, Phys. Rev. 
E 71, 011105 (2005). 

R. D. Kamien and A. J. Liu, Phys. Rev. Lett. 99, 155501 
(2007). 

We use the canonical notation for the dynamical transi- 
tion ^d, but the reader should not confuse this d sub- 
script with dimensionality d. 

R. Mari and J. Kurchan, arXiv: 1104.3420 (2011). 
G. Brambilla et al. , Phys. Rev. Lett. 102, 085703 (2009) 



[30] Within RFOT, the "maximally random jammed" pack- 
ing discussed in Ref. [ITI ] should be very close to ipth UM- 

[31] F. Krzakala and J. Kurchan, Phys. Rev. E 76, 021122 
(2007); S. Torquato and F. H. StiUinger, Rev. Mod. Phys. 
82, 2633 (2010). 

[32] A. Andreanov, G. Biroli, and J. -P. Bouchaud, Europhys. 
Lett. 88, 16001 (2009). 



J 



[33] H. Jacquin and F. van Wijland, Phys. Rev. Lett. 106, 
210602 (2011); Y. Elskens and H. L. Frisch, Phys. Rev. 
A 37, 4351 (1988) . 

[34] M. Bishop, N. Clisby, and P. A. Whitlock, J. Chem. 
Phys. 128, 034506 (2008). 



APPENDIX 



d 


N 


Af 


A^/ 


1/pmin 


/th 


2Vd 


2Vth 


2'^^K 


2 (fiGCP 


2 V5MCT 








Ref. [24] 










(RT+CS) (RT-FCS) (MCT-hPY) 
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0.5 
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0.03 


0.11(2) 


4.57(6) 


5.21(1) 
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5.47 


4.13 
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0.114 


-0.051 


0.02 


0.12(1) 


6.41(6) 


7.48(1) 


6.91 


7.79 
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-1.086 


-1.625 
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-4.910 
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-35.46 


-41.10 


0.02 


0.137(4) 
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25.51(7) 
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27.4 


22.2 
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20000 


-62.50 


-73.83 


0.02 


0.139(4) 


23.1(2) 


30.37(8) 


27.2 


32.7 


27.0 
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50000 


-107.0 


-129.6 


0.015 


0.141(1) 


26.9(2) 


35.64(9) 


31.7 


38.3 


32.4 


12 
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-176.0 
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44.2 
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TABLE IL Numerical values of the RFOT critical densities and of the simulation parameters. The simulation results are 
extrapolated at 7 = according to the procedure discussed in the text. 



Error analysis - From the simulation results, the error on A^' is found to be very small and, for all practical 
purposes, the quantity is here considered exact. While the error on (^th and /th are estimated by taking three times 
the error on the linear regression coefficients, for (/jj, a standard error propagation is conducted starting from its 
definition. Indeed (pd is the intersection of the fluid EOS and the threshold glass EOS defined by the extrapolation of 
the free volume equation to 7 = 0: 



Defining q{ip) = ^/ptiuidiv), we then get that 
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dipthq{ipd) 



l + d{l- fth)vthq'{p'd) 



Sfi 
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Finite-size effects - The finite-size arguments that we presented in the main text arc tested by comparing d = 8 
systems with A^ = 1.4K, 2.5K, 8K and 32K. Although the number of particles increases by nearly 50% per dimension 
over this sequence, the compaction curves are almost indistinguishable except for a slight finite size effect for the 
largest 7 = 0.1. The results of our data analysis coincide within numerical errors, estimated as discussed in the 
previous paragraph (Table HI]) . 



